// always include the config file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

// C++ includes
#include<math.h>
#include<iostream>
#include<vector>
#include<map>
#include<string>
#include<sys/stat.h>

// dune-common includes
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/common/timer.hh>
#include<dune/common/exceptions.hh>
#include<dune/common/fvector.hh>
#include<dune/common/typetraits.hh>
#include <dune/common/float_cmp.hh>

// dune-grid includes
#include<dune/grid/io/file/vtk.hh>
#include<dune/grid/io/file/vtk/vtkwriter.hh>
#include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include<dune/grid/io/file/gmshreader.hh>
#include<dune/grid/onedgrid.hh>
#include<dune/grid/yaspgrid.hh>
#include<dune/grid/utility/structuredgridfactory.hh>
#if HAVE_UG
#include<dune/grid/uggrid.hh>
#endif
#if HAVE_DUNE_ALUGRID
#include<dune/alugrid/grid.hh>
#include<dune/alugrid/dgf.hh>
#include<dune/grid/io/file/dgfparser/dgfparser.hh>
#endif

// dune-geometry includes
#include<dune/geometry/referenceelements.hh>
#include<dune/geometry/quadraturerules.hh>

// dune-istl included by pdelab
// dune-pdelab includes
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/common/vtkexport.hh>
#include<dune/pdelab/common/instationaryfilenamehelper.hh>

#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include<dune/pdelab/gridfunctionspace/interpolate.hh>
#include<dune/pdelab/gridfunctionspace/vtk.hh>
#include<dune/pdelab/gridfunctionspace/subspace.hh>
#include<dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>

#include<dune/pdelab/finiteelementmap/pkfem.hh>
#include<dune/pdelab/finiteelementmap/qkfem.hh>
#include<dune/pdelab/finiteelementmap/qkdg.hh>
#include<dune/pdelab/constraints/common/constraints.hh>
#include<dune/pdelab/constraints/common/constraintsparameters.hh>
#include<dune/pdelab/constraints/conforming.hh>

#include<dune/pdelab/gridoperator/gridoperator.hh>
#include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
#include<dune/pdelab/localoperator/variablefactories.hh>

#include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/function/callableadapter.hh>
#include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/stationary/linearproblem.hh>
#include<dune/pdelab/solver/newton.hh>

#include<dune/istl/bvector.hh>
#include<dune/istl/operators.hh>
#include<dune/istl/solvers.hh>
#include<dune/istl/preconditioners.hh>
#include<dune/istl/io.hh>
#include<dune/istl/superlu.hh>

//===============================================================
// Include your hyperbolic model and problem to solve
//===============================================================

#include"hyperbolicdg.hh"
#include"shallowwater.hh" //model
#include"shallowwaterproblem.hh"//problem
#include"numericalflux.hh"
#include"driver.hh"

//===============================================================
// Main program with grid setup
//===============================================================
int main(int argc, char** argv)
{
  try{
    // Maybe initialize Mpi
    Dune::MPIHelper&
      helper = Dune::MPIHelper::instance(argc, argv);
    if(Dune::MPIHelper::isFake)
      std::cout<< "This is a sequential program." << std::endl;
    else
      std::cout << "Parallel code run on "
                << helper.size() << " process(es)" << std::endl;

    // open ini file
    Dune::ParameterTree ptree;
    Dune::ParameterTreeParser ptreeparser;
    ptreeparser.readINITree("step00-shallowwater-mpi.ini",ptree);
    ptreeparser.readOptions(argc,argv,ptree);

    // read ini file
    const int dim = ptree.get<int>("grid.dim");
    const int refinement = ptree.get<int>("grid.refinement");
    std::string gridmanager = ptree.get("grid.manager","yasp");    
    const int degree = ptree.get<int>("fem.degree");

    // in 1d use OneDGrid: only sequence version, no parallel.
    // // YaspGrid section
    // if (dim==1 && gridmanager=="yasp")
    //   {
    //     const int dim=1;
    //     typedef Dune::YaspGrid<dim> Grid;
    //     typedef Grid::ctype DF;
    //     Dune::FieldVector<DF,dim> L;
    //     L[0] = ptree.get("grid.structured.LX",(double)1.0);
    //     std::array<int,dim> N;
    //     N[0] = ptree.get("grid.structured.NX",(int)200);
    //     std::bitset<dim> B(false);
    //     int overlap=1;
    //     std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>
    //       (new Grid(L,N,B,overlap,Dune::MPIHelper::getCollectiveCommunication()));
    //     gridp->refineOptions(false); // keep overlap in cells
    //     gridp->globalRefine(refinement);
    //     typedef Grid::LeafGridView GV;
    //     GV gv=gridp->leafGridView();
    //     // create problem (setting)
    //     using PROBLEM = Problem<GV,GV::Grid::ctype>;
    //     PROBLEM problem;
    //     // create model on a given setting
    //     using MODEL = Model<PROBLEM>;
    //     MODEL model(problem);
    //     // create numerical flux
    //     using NUMFLUX = LLFflux<MODEL>;
    //     NUMFLUX numflux(model);        
    //     if (degree==0) {
    //       typedef Dune::PDELab::QkDGLocalFiniteElementMap
    //         <GV::Grid::ctype,double,0,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
    //       FEM fem;
    //       driver(gv,fem,numflux,ptree);
    //     }
    //     if (degree==1) {
    //       typedef Dune::PDELab::QkDGLocalFiniteElementMap
    //         <GV::Grid::ctype,double,1,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
    //       FEM fem;
    //       driver(gv,fem,numflux,ptree);
    //     }        
    //     if (degree==2) {
    //       typedef Dune::PDELab::QkDGLocalFiniteElementMap
    //         <GV::Grid::ctype,double,2,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
    //       FEM fem;
    //       driver(gv,fem,numflux,ptree);
    //     }
    //   }
    // YaspGrid section
    if (dim==2 && gridmanager=="yasp")
      {
        const int dim=2;
        typedef Dune::YaspGrid<dim> Grid;
        typedef Grid::ctype DF;
        Dune::FieldVector<DF,dim>maxCoord;
        maxCoord[0] = ptree.get("gridYasp.MAXX",(double)1.0);
        maxCoord[1] = ptree.get("gridYasp.MAXY",(double)1.0);
        std::array<int,dim> elementN;
        elementN[0] = ptree.get("gridYasp.NX",(int)1);
        elementN[1] = ptree.get("gridYasp.NY",(int)1);

        std::bitset<dim> Boun(false);
        int overlap=1;
        std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>
          (new Grid(maxCoord,elementN,Boun,overlap,helper.getCommunicator()));
		    // Dune::MPIHelper::getCollectiveCommunication()));

        gridp->refineOptions(false); // keep overlap in cells
        gridp->globalRefine(refinement);
        typedef Grid::LeafGridView GV;
        GV gv=gridp->leafGridView();
        // create problem (setting)
        using PROBLEM = Problem<GV,GV::Grid::ctype>;
        PROBLEM problem;
        // create model on a given setting
        using MODEL = Model<PROBLEM>;
        MODEL model(problem);
        // create numerical flux
        using NUMFLUX = LLFflux<MODEL>;
        NUMFLUX numflux(model);        
        if (degree==0) {
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,0,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          driver(gv,fem,numflux,ptree);
        }
        if (degree==1) {
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,1,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          driver(gv,fem,numflux,ptree);
        }        
        if (degree==2) {
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,2,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          driver(gv,fem,numflux,ptree);
        }
      }
    // UG section    
    if (dim==2 && gridmanager=="ug")
      {
# if HAVE_UG
        const int dim=2;	
        typedef Dune::UGGrid<dim> Grid;
        std::string filename = ptree.get("grid.twod.filename",
                                         "unitsquare.msh");
        Dune::GridFactory<Grid> factory;	
	// std::cout << "error location 0" << std::endl;
        Dune::GmshReader<Grid>::read(factory,filename,true,true);
	// std::cout << "error location 1" << std::endl;
        std::shared_ptr<Grid> gridp(factory.createGrid());
	std::cout << "error location 2" << std::endl;	
        Dune::Timer timer;
        gridp->globalRefine(refinement);
        gridp->loadBalance();
        std::cout << "Time for mesh refinement " << timer.elapsed()
                  << " seconds" << std::endl;
        typedef Grid::LeafGridView GV;
        typedef Grid::ctype DF;
        GV gv=gridp->leafGridView();
        // create problem (setting)
        using PROBLEM = Problem<GV,GV::Grid::ctype>;
        PROBLEM problem;
        // create model on a given setting
        using MODEL = Model<PROBLEM>;
        MODEL model(problem);
        // create numerical flux
        using NUMFLUX = LLFflux<MODEL>;
        NUMFLUX numflux(model);
	std::cout << "error location 3" << std::endl;		
        if (degree==0) {
	  typedef Dune::PDELab::QkDGLocalFiniteElementMap
	    <GV::Grid::ctype,double,0,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
	  FEM fem;
          driver(gv,fem,numflux,ptree);
        }
        if (degree==1) {
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,1,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
	std::cout << "error location 4" << std::endl;			  
          driver(gv,fem,numflux,ptree);
        }        
        if (degree==2) {
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,2,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          driver(gv,fem,numflux,ptree);
        }
	std::cout << "error location 5" << std::endl;		
#else
        std::cout << "UG was not found during installation" << std::endl;
#endif        
      }
    // YASP section
    if (dim==3 && gridmanager=="yasp")
      {
        const int dim=3;
        typedef Dune::YaspGrid<dim> Grid;
        typedef Grid::ctype DF;
        Dune::FieldVector<DF,dim> L;
        L[0] = ptree.get("grid.structured.LX",(double)1.0);
        L[1] = ptree.get("grid.structured.LY",(double)1.0);
        L[2] = ptree.get("grid.structured.LZ",(double)1.0);
        std::array<int,dim> N;
        N[0] = ptree.get("grid.structured.NX",(int)16);
        N[1] = ptree.get("grid.structured.NY",(int)16);
        N[2] = ptree.get("grid.structured.NZ",(int)16);
    //     std::bitset<dim> B(false);
    //     int overlap=1;
    //     std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>(new Grid(L,N,B,overlap,Dune::MPIHelper::getCollectiveCommunication()));
    //     gridp->refineOptions(false); // keep overlap in cells
    //     gridp->globalRefine(refinement);
    //     typedef Grid::LeafGridView GV;
    //     GV gv=gridp->leafGridView();
    //     // create problem (setting)
    //     using PROBLEM = Problem<GV,GV::Grid::ctype>;
    //     PROBLEM problem;
    //     // create model on a given setting
    //     using MODEL = Model<PROBLEM>;
    //     MODEL model(problem);
    //     // create numerical flux
    //     using NUMFLUX = LLFflux<MODEL>;
    //     NUMFLUX numflux(model);        
    //     if (degree==0) {
    //       typedef Dune::PDELab::QkDGLocalFiniteElementMap
    //         <GV::Grid::ctype,double,0,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
    //       FEM fem;
    //       driver(gv,fem,numflux,ptree);
    //     }
    //     if (degree==1) {
    //       typedef Dune::PDELab::QkDGLocalFiniteElementMap
    //         <GV::Grid::ctype,double,1,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
    //       FEM fem;
    //       driver(gv,fem,numflux,ptree);
    //     }        
    //     if (degree==2) {
    //       typedef Dune::PDELab::QkDGLocalFiniteElementMap
    //         <GV::Grid::ctype,double,2,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
    //       FEM fem;
    //       driver(gv,fem,numflux,ptree);
    //     }
      }
//     // UG section
//     if (dim==3 && gridmanager=="ug")
//       {
// # if HAVE_UG
//         typedef Dune::UGGrid<3> Grid;
//         std::string filename = ptree.get("grid.threed.filename",
//                                          "unitsquare.msh");
//         Dune::GridFactory<Grid> factory;
//         Dune::GmshReader<Grid>::read(factory,filename,true,true);
//         std::shared_ptr<Grid> gridp(factory.createGrid());
//         Dune::Timer timer;
//         gridp->globalRefine(refinement);
//         gridp->loadBalance();        
//         std::cout << "Time for mesh refinement " << timer.elapsed()
//                   << " seconds" << std::endl;
//         typedef Grid::LeafGridView GV;
//         typedef Grid::ctype DF;
//         GV gv=gridp->leafGridView();
//         // create problem (setting)
//         using PROBLEM = Problem<GV,GV::Grid::ctype>;
//         PROBLEM problem;
//         // create model on a given setting
//         using MODEL = Model<PROBLEM>;
//         MODEL model(problem);
//         // create numerical flux
//         using NUMFLUX = LLFflux<MODEL>;
//         NUMFLUX numflux(model);
//         if (degree==0) {
//           typedef Dune::PDELab::QkDGLocalFiniteElementMap
//             <GV::Grid::ctype,double,0,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
//           FEM fem;
//           driver(gv,fem,numflux,ptree);
//         }
//         if (degree==1) {
//           typedef Dune::PDELab::QkDGLocalFiniteElementMap
//             <GV::Grid::ctype,double,1,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
//           FEM fem;
//           driver(gv,fem,numflux,ptree);
//         }        
//         if (degree==2) {
//           typedef Dune::PDELab::QkDGLocalFiniteElementMap
//             <GV::Grid::ctype,double,2,dim,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
//           FEM fem;
//           driver(gv,fem,numflux,ptree);
//         }
// #else
//         std::cout << "UG was not found during installation" << std::endl;
// #endif        
//       }
  }
  catch (Dune::Exception &e){
    std::cerr << "Dune reported error: " << e << std::endl;
    return 1;
  }
  catch (...){
    std::cerr << "Unknown exception thrown!" << std::endl;
    return 1;
  }
}
